# Acting under Chapter 7: Analysis #
# Author: Johannes Scherzinger #
# Date: 08.02.2021 #

# Packages
library(stargazer)
library(survival)

# To run the models below, all independent variables (and the DV) have to be created from Script
# using "ActingChap7_processing.R" (08.02.21)



# run simple logit model on rhetoric and authorization
logit.model <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$gen.pres) + as.numeric(UNSC$r2p.pres) + as.numeric(UNSC$crisis.pres) + as.numeric(UNSC$hr.pres) + as.numeric(UNSC$terror.pres) + as.numeric(UNSC$dem.pres) + as.numeric(UNSC$reg.pres), family = binomial(link = 'logit'), data = UNSC)
summary(logit.model)

# convert to odds ratios
exp(coef(logit.model))

# get OR and CIs in one table
exp(cbind(OR = coef(logit.model), confint(logit.model)))


# we have reason to believe that time-invariant change may affect the outcome
# therefore we can use a within-change estimator to tackle this sort of unobserved bias 
# in this instance on country level
# note that "R" has no standard package for fixed effects models
#year dummies are possible but clogg output
#you can use the clogit function from the survival packages which the authors have recoded to do the same

# run logit country fixed effects
logit.fixed.country <- clogit(as.numeric(UNSC$auth) ~ as.numeric(UNSC$gen.pres) + as.numeric(UNSC$r2p.pres) + as.numeric(UNSC$crisis.pres) + as.numeric(UNSC$hr.pres) + as.numeric(UNSC$terror.pres) + as.numeric(UNSC$dem.pres) + as.numeric(UNSC$reg.pres) + strata(UNSC$country), data = UNSC)
summary(logit.fixed.country) # note that this might take a while depending on your machine



# estimate a model with a p5 control variable
# run logit controlling for p5
logit.p5 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$gen.pres) + as.numeric(UNSC$r2p.pres) + as.numeric(UNSC$crisis.pres) + as.numeric(UNSC$hr.pres) + as.numeric(UNSC$terror.pres) + as.numeric(UNSC$dem.pres) + as.numeric(UNSC$reg.pres) + as.factor(UNSC$P5), family = binomial(link = 'logit'), data = UNSC)
summary(logit.p5)

# convert to odds ratios
exp(coef(logit.p5))



# run logit conditional p5, year-fixed-effects
logit.p5.fixed <- clogit(as.numeric(UNSC$auth) ~ as.numeric(UNSC$gen.pres) + as.numeric(UNSC$r2p.pres) + as.numeric(UNSC$crisis.pres) + as.numeric(UNSC$hr.pres) + as.numeric(UNSC$terror.pres) + as.numeric(UNSC$dem.pres) + as.numeric(UNSC$reg.pres) + strata(year.x) + as.factor(UNSC$P5), data = UNSC)
summary(logit.p5.fixed)



# export  models with stargazer (Hvalec 2018)
stargazer(logit.model, logit.fixed.country, logit.p5, logit.p5.fixed, type = "html", style = "APSR", out = "testmodel_final.htm")







##### APPENDIX #####


# estimate a model with time polynomials 
logit.poly <- glm(as.numeric(UNSC$auth) ~  as.numeric(UNSC$gen.pres) + as.numeric(UNSC$r2p.pres) + as.numeric(UNSC$crisis.pres) + as.numeric(UNSC$hr.pres) + as.numeric(UNSC$terror.pres) + as.numeric(UNSC$dem.pres) + as.numeric(UNSC$reg.pres) + UNSC$year.x + I(UNSC$year.x^2) + I(UNSC$year.x^3), data = UNSC)
summary(logit.poly)

# convert to odds ratios
exp(coef(logit.poly))

# estimate additional bivariate regressions #
bv1 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$gen.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv1)

bv2 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$r2p.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv2)

bv3 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$crisis.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv3)

bv4 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$hr.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv4)

bv5 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$terror.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv5)

bv6 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$dem.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv6)

bv7 <- glm(as.numeric(UNSC$auth) ~ as.numeric(UNSC$reg.pres), family = binomial(link = 'logit'), data = UNSC)
summary(bv7)


# export models with stargazer
stargazer(logistic.model, logit.fixed, logit.poly, type = "html", style = "APSR", out = "testmodel.htm")

# export even more models
stargazer(logistic.model, logit.fixed, logistic.p5, logistic.e10, logit.poly,type = "html", style = "APSR", out = "testmodel3.htm")


# export models in Appendix
stargazer(bv1,bv2,bv3,bv4,bv5,bv6,bv7, logit.poly, type = "html", style = "APSR", out = "appendix_model.htm")



############ Replicate Analysis using RESOLUTION DEBATES #######################
# Aggregate speeches to RESOLUTION DEBATES , deflate DV

UNSC_2 <- UNSC %>% 
  select(country, everything()) %>% 
  group_by(resnum,year.x) %>% 
  summarize(mean.gen = mean(gen.pres),
            mean.hr = mean(hr.pres),
            mean.cri = mean(crisis.pres),
            mean.reg = mean(reg.pres),
            mean.dem = mean(dem.pres),
            mean.terror = mean(terror.pres),
            mean.r2p = mean(r2p.pres),
            mean.auth = mean(auth))


# run again
logit.2 <- glm(mean.auth ~ mean.gen + mean.hr + mean.cri + mean.r2p +
                 mean.reg + mean.dem + mean.terror, family = binomial, data = UNSC_2)
summary(logit.2)


# run controlling for year-fixed-effects
# using plm package
library(plm)
logit.3 <- plm(mean.auth ~ mean.gen + mean.hr + mean.terror+
                 mean.cri + mean.reg+ mean.dem + mean.r2p, index = "year.x", model = "within", data = UNSC_2)

summary(logit.3)
exp(cbind(OR = coef(logit.fixed.country)))


# finally, use lagged dependent variable to reduce possible bias in anticipated resolution outcome
# lag outcome use Hmisc Package
library(Hmisc)

# Lag DV
UNSC_2$l.mean.auth <- Lag(UNSC_2$mean.auth)

# now estimate linear probability model | in "R" its simply lm
linear.prob <-lm(mean.auth ~ mean.gen + mean.hr + mean.cri + mean.r2p+ mean.reg +
                   mean.dem + mean.terror + l.mean.auth, data = UNSC_2)
summary(linear.prob)


# Export Again using stargazer
stargazer(logit.2, logit.3, linear.prob, type = "html", style = "APSR", out = "appendix_debates_models.htm")

